library(scran)
library(dplyr)
library(tidyr)
library(ggplot2)
library(scuttle)
library(pheatmap)
library(patchwork)
library(stringr)
library(Seurat)
library(ggbeeswarm)
library(EnhancedVolcano)
library(SingleCellExperiment)
library(gridExtra)DE_LECs
Run Differential expression analysis within LEC subtypes
Standard DE analysis using pseudobulk DE methods (note this differs from the FGCZ default analysis). We want to identify differentially expressed genes within LEC cluster between all groups of skin, YUMMER and YUMM.
Preamble
Code
# volcano plot function as defined in https://github.com/HelenaLC/TLS-Silina/blob/main/code/geo-02-differential.qmd
.volcano <- \(df, title, fdr = 0.05, lfc = 1, select_lab = NULL) {
EnhancedVolcano(df,
x = "logFC", y = "PValue",
FCcutoff = lfc, pCutoff = fdr,
selectLab = select_lab,
pointSize = 1.7, raster = TRUE,
title = title, subtitle = NULL,
lab = rownames(df), labSize = 4,
drawConnectors = TRUE, widthConnectors = 0.5) +
guides(col = guide_legend(override.aes = list(alpha = 1, size = 5))) +
theme_bw(9) + theme(
aspect.ratio = 1,
legend.title = element_blank(),
panel.grid.minor = element_blank())
}
.heatmap <- \(mtx, de_genes, cd, title, fdr = 0.05, n_lfc = 20) {
top <- de_genes %>%
filter(PValue < fdr) |>
slice_max(abs(logFC), n = n_lfc)
mtx_sub <- log1p(mtx[rownames(top),])
if (length(rownames(top)) < 2){
return(print("No de genes"))
}else{
hm <- pheatmap(mtx_sub,
main = title, fontsize = 6,
col = rev(hcl.colors(51, "RdBu")),
scale = "row", show_colnames = FALSE, annotation_col = cd)
hm
}
}Data object
seurat<- readRDS(file.path("..", "..","..","data", "scData_LEC_tumor_skin.rds"))
# correct condition assignment!!
seurat$cond <- seurat[[]] |>
mutate(
cond = case_when(
str_detect(Sample, "YUMM[0-9]") ~ "YUMM",
str_detect(Sample, "YUMMER") ~ "YUMMER",
str_detect(Sample, "Skin") ~ "skin"
)
) |> select(cond)
# check assignment
table(seurat$Sample, seurat$cond)
skin YUMM YUMMER
SkinLECs_Leukocytes1 469 0 0
SkinLECs_Leukocytes2 669 0 0
TumorYUMM1_1A 0 17 0
TumorYUMM1_1B 0 18 0
TumorYUMM2_1A 0 57 0
TumorYUMM2_1B 0 49 0
TumorYUMM5_2A 0 44 0
TumorYUMM5_2B 0 44 0
TumorYUMM6_2A 0 58 0
TumorYUMM6_2B 0 72 0
TumorYUMMER3_1A 0 0 66
TumorYUMMER3_1B 0 0 44
TumorYUMMER4_1A 0 0 81
TumorYUMMER4_1B 0 0 68
TumorYUMMER7_2A 0 0 86
TumorYUMMER7_2B 0 0 101
TumorYUMMER8_2A 0 0 160
TumorYUMMER8_2B 0 0 135seurat$sample_pooled <- seurat$Sample |>
forcats::fct_collapse("TumorYUMM1" = c("TumorYUMM1_1A", "TumorYUMM1_1B"),
"TumorYUMM2" = c("TumorYUMM2_1A", "TumorYUMM2_1B"),
"TumorYUMM5" = c("TumorYUMM5_2A", "TumorYUMM5_2B"),
"TumorYUMM6" = c("TumorYUMM6_2A", "TumorYUMM6_2B"))
table(seurat$cond)
skin YUMM YUMMER
1138 359 741 # switch to SingleCellExperiment object
sce <- as.SingleCellExperiment(seurat)
sce$cond <- as.factor(sce$cond)DE analysis
cond_de <- function(cond1, cond2){
# subset sce
sce_sub <- sce[,sce$cond %in% c(cond1, cond2)]
sce_sub$cond <- droplevels(sce_sub$cond)
# creating pseudobulks
summed <- aggregateAcrossCells(sce_sub,
id=colData(sce_sub)[,c("ident", "Sample")])
# filter min cells
summed.filt <- summed[,summed$ncells >= 5]
print(table(summed.filt$cond, summed.filt$ident))
# de
design <- model.matrix(~ summed.filt$cond)
de.results <- pseudoBulkDGE(summed.filt,
label=summed.filt$ident,
design=~factor(cond),
coef=ncol(design),
condition=summed.filt$cond
)
}
yumm_skin <- cond_de("skin", "YUMM")
0 1 2 3 4 5 6 7
skin 2 2 2 2 2 0 2 1
YUMM 6 4 5 3 3 6 1 5yummer_yumm <- cond_de("YUMMER", "YUMM")
0 1 2 3 4 5 6 7
YUMM 6 4 5 3 3 6 1 5
YUMMER 8 5 8 8 6 7 5 2yummer_skin <- cond_de("skin","YUMMER")
0 1 2 3 4 5 6 7
skin 2 2 2 2 2 0 2 1
YUMMER 8 5 8 8 6 7 5 2# refactor cond
sce$cond <- factor(sce$cond, levels = c("YUMMER", "YUMM", "skin"))
skin_yumm <- cond_de("skin", "YUMM")
0 1 2 3 4 5 6 7
YUMM 6 4 5 3 3 6 1 5
skin 2 2 2 2 2 0 2 1#DE skin_yumm
is.de <- decideTestsPerLabel(skin_yumm, threshold=0.05)
summarizeTestsPerLabel(is.de) -1 0 1 NA
0 56 1216 79 11727
1 11 538 13 12516
2 26 773 64 12215
3 32 444 52 12550
4 8 623 8 12439
6 42 1663 59 11314
7 131 749 44 12154#DE yumm_skin
is.de <- decideTestsPerLabel(yumm_skin, threshold=0.05)
summarizeTestsPerLabel(is.de) -1 0 1 NA
0 79 1216 56 11727
1 13 538 11 12516
2 64 773 26 12215
3 52 444 32 12550
4 8 623 8 12439
6 59 1663 42 11314
7 44 749 131 12154#DE yummer_yumm
is.de <- decideTestsPerLabel(yummer_yumm, threshold=0.05)
summarizeTestsPerLabel(is.de) -1 0 1 NA
0 1 1155 15 11907
1 1 453 2 12622
2 6 999 11 12062
3 0 1985 13 11080
4 1 722 4 12351
5 3 529 10 12536
6 9 1403 1 11665
7 0 670 4 12404#DE yummer_skin
is.de <- decideTestsPerLabel(yummer_skin, threshold=0.05)
summarizeTestsPerLabel(is.de) -1 0 1 NA
0 193 2059 168 10658
1 29 1262 26 11761
2 216 1755 152 10955
3 284 2497 297 10000
4 10 992 18 12058
6 114 1566 93 11305
7 44 456 94 12484#write as output
write_de_res <- function(de_obj, res_nam){
res_list <- lapply(names(de_obj), function(clust, nam = res_nam){
cur.res <- de_obj[[clust]]
res <- cur.res[order(cur.res$PValue),] |>
as.data.frame() |>
filter(!is.na(PValue))
## --- adapt path until DE/ to local output path -----##
write.csv(res, paste0("../../../out/DE/LEC_", nam, "_", clust, ".csv"))
res
})
}
de_yummer_yumm <- write_de_res(yummer_yumm, res_nam = "yummer_yumm")
de_yumm_skin <- write_de_res(yumm_skin, res_nam = "yumm_skin")
de_yummer_skin <- write_de_res(yummer_skin, res_nam = "yummer_skin")Plot results
yumm vs skin
plot_list <- lapply(names(yumm_skin), \(.) {
df = yumm_skin[[.]]
.volcano(df = df, title = ., fdr = 0.05, lfc = 1)
}
)
wrap_plots(plot_list, nrow = 4) +
plot_layout(guides = "collect") &
theme(legend.position = "top")
skin vs yumm
plot_list <- lapply(names(skin_yumm), \(.) {
df = skin_yumm[[.]]
.volcano(df = df, title = ., fdr = 0.05, lfc = 1)
}
)
wrap_plots(plot_list, nrow = 4) +
plot_layout(guides = "collect") &
theme(legend.position = "top")
yummer vs yumm
plot_list <- lapply(names(yummer_yumm), \(.) {
df = yummer_yumm[[.]]
.volcano(df = df, title = ., fdr = 0.05, lfc = 1)
}
)
wrap_plots(plot_list, nrow = 4) +
plot_layout(guides = "collect") &
theme(legend.position = "top")
yummer vs skin
plot_list <- lapply(names(yummer_skin), \(.) {
df = yummer_skin[[.]]
.volcano(df = df, title = ., fdr = 0.05, lfc = 1, select_lab = c("Lyve1"))
.volcano(df = df, title = ., fdr = 0.05, lfc = 1)
}
)
wrap_plots(plot_list, nrow = 4) +
plot_layout(guides = "collect") &
theme(legend.position = "top")
Heatmap plots
yumm vs skin
for (. in names(yumm_skin)) {
cat("####", ., "\n")
sub <- subset(x = seurat, idents = .)
sub <- subset(x = sub, subset = cond %in% c("YUMM", "skin"))
mtx <- GetAssayData(object = sub, slot = 'data')
top <- as.data.frame(yumm_skin[[.]]) %>%
filter(PValue < 0.05) |>
slice_max(abs(logFC), n = 40)
mtx_sub <- log1p(mtx[rownames(top),])
cd <- data.frame(sub[[]] |> select('cond'))
hm <- pheatmap(mtx_sub,
main = ., fontsize = 6,
col = rev(hcl.colors(51, "RdBu")),
scale = "row",
show_colnames = FALSE,
cluster_cols = FALSE,
annotation_col = cd)
print(hm); cat("\n\n")
}#### 0 
#### 1 
#### 2 
#### 3 
#### 4 
#### 6 
#### 7 
yummer vs yumm
for (. in names(yummer_yumm)) {
cat("####", ., "\n")
sub <- subset(x = seurat, idents = .)
sub <- subset(x = sub, subset = cond %in% c("YUMMER", "YUMM"))
mtx <- GetAssayData(object = sub, slot = 'data')
top <- as.data.frame(yummer_yumm[[.]]) %>%
filter(PValue < 0.05) |>
slice_max(abs(logFC), n = 40)
mtx_sub <- log1p(mtx[rownames(top),])
cd <- data.frame(sub[[]] |> select('cond'))
hm <- pheatmap(mtx_sub,
main = ., fontsize = 6,
col = rev(hcl.colors(51, "RdBu")),
scale = "row",
show_colnames = FALSE,
cluster_cols = FALSE,
annotation_col = cd)
print(hm); cat("\n\n")
}#### 0 
#### 1 
#### 2 
#### 3 
#### 4 
#### 5 
#### 6 
#### 7 
yummer vs skin
for (. in names(yummer_skin)) {
cat("####", ., "\n")
sub <- subset(x = seurat, idents = .)
sub <- subset(x = sub, subset = cond %in% c("YUMMER", "skin"))
mtx <- GetAssayData(object = sub, slot = 'data')
top <- as.data.frame(yummer_skin[[.]]) %>%
filter(PValue < 0.05) |>
slice_max(abs(logFC), n = 40)
mtx_sub <- log1p(mtx[rownames(top),])
cd <- data.frame(sub[[]] |> select('cond'))
hm <- pheatmap(mtx_sub,
main = ., fontsize = 6,
col = rev(hcl.colors(51, "RdBu")),
scale = "row",
show_colnames = FALSE,
cluster_cols = FALSE,
annotation_col = cd)
print(hm); cat("\n\n")
}#### 0 
#### 1 
#### 2 
#### 3 
#### 4 
#### 6 
#### 7 